Pressure-energy correlations in liquids. III. Statistical mechanics and thermodynamics of liquids 

with hidden scale invariance 



Thomas B. Schr0der,''[2] Nicholas P. Bailey,' [tjuif R. Pedersen,' [t]Nicoletta Gnan,''[§]and Jeppe C. Dyre'-j^ 

' DNRF Center "Glass and Time", IMFUFA, Dept. of Sciences, 
Roskilde University, P.O. Box 260, DK-4000 Roskilde, Denmark 
(Dated: August 21, 2009) 

In this third paper of the series, which started with [N. P. Bailey et al, J. Chem. Phys. 129, 184507 and 
184508 (2008)], we continue the development of the theoretical understanding of strongly correlating liquids 
- those whose instantaneous potential energy and virial are strongly correlated in their thermal equilibrium 
fluctuations at constant volume. The existence of such liquids was detailed in previous work which identified 
them, based on computer simulations, as a large class of liquids, including van der Waals liquids but not, 
e.g., hydrogen-bonded liquids. We here discuss the following: (1) The scaling properties of inverse power- 
law and extended inverse power-law potentials (the latter include a linear term which "hides" the approximate 
scale invariance); (2) results from computer simulations of molecular models concerning out-of-equilibrium 
conditions; (3) ensemble dependence of the virial / potential energy correlation coefficient; (4) connection to 
the Griineisen parameter; (5) interpretation of strong correlations in terms of the energy-bond formalism. 



I. INTRODUCTION 



In a series of papers published last yeapDSEEEl introduced the concept of strongly correlating liquids and demonstrated by 
computer simulations that this includes a large class of model liquids. Specifically, the fluctuations which are in many cases 
strongly coiTelated are those of the configurational parts of pressure and energy, i.e., the parts in addition to the ideal gas terms, 
coming from the interatomic forces. Recall that for any microscopic state, energy E and pressure p have contributions both from 
particle momenta and positions: 

E = K{pi, Pat) + C/(ri, . . . , rjv) 

p = iVfcsr(pi, . . . , pn)/V + W{ri, . . .,rN)/V . (1) 

Here K and U are the kinetic and potential energies, respectively, and T{pi, . . . , pjv) is the "kinetic temperature", proportional 
to the kinetic energy per particle.- The configurational contribution to pressure is the virial W, which is defined^ by 

(2) 

i 

For a liquid with pair interactions, if v{r) is the pair potential and Vij is the distance between particles i and j, we have 



Upidr ^^v{rij) (3) 

i<j 



Strong W, U correlation, if present at all, is observed under conditions of fixed volume, as illustrated in Fig.[TJa). The degree of 
correlation is quantified by the standard correlation coefficient R, defined^ by 



(AWAU) 

Here and henceforth, unless otherwise specified, angle brackets () denote thermal NVT ensemble averages; A denotes deviation 
from the average value of the quantity in qu estion. We call liquids with R > 0.9 strongly correlating. Another characteristic 
quantity is the "slope" 7, which we define^as the ratio of standard deviations: 
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In the limit of perfect correlation (i? ^ 1) 7 becomes equal to the standard linear-regression slope for W as a function of U at 
fixed volume. 

In Paper P of this series it was shown that strongly correlating liquids are typically those with van der Waals type attraction 
and steep repulsion, which in simulations are often modelled by combinations of one or more Lennard-Jones type potentials. 
Typical slope values for the latter are of order 6, depending slightly on state point (in the limit of very high density or temperature 
the slope converges slowly to 4). Experimental data for Argon were analyzed and shown to be consistent with strong correlations 
(R > 0.95) in the region of the phase diagram where quantum effects are not important. 

It is worth noting that the class of strongly correlating liquids does not simply correspond to radially symmetric pair potentials. 
Firstly, two metallic systems with many-body potentials were found to be strongly correlating; it is probably true for metallic 
systems in general, although this needs to be confirmed. Also many molecular liquids are strongly correlating. In fact, any 
potential with an inverse power-law dependence on distances (not necessarily based on pair interactions) is perfectly correlating. 
Secondly, there exist radially symmetric pair potentials which are not strongly correlating, for example the Dzugutov system. "^El 
One reason for strong correlation not to hold in some molecular systems is the presence of Coulombic terms in the potential. 
By themselves these would give strong correlation, but their combination with Lennard-Jones forces typically leads to weak 
correlation. This was detailed in Paper I, which presented results from simulations of 13 different model liquids. In our present 
understanding based on these simulations, liquids with two length scales in their potentials are rarely strongly correlating. 

Paper 10 in this series analyzed the case of the standard single-component Lennard-Jones liquid in detail. Building on the 
fact that inverse power-law potentials cx r~" are perfectly correlating, the results of this analysis can be summarized as follows: 
(1) Almost all of the fluctuations in W and U come from interparticle separations in the region of the first peak of the radial 
distribution function g{r); (2) in this region the Lennard-Jones potential is approximated very well by the sum of an inverse 
power law with exponent n ^18 and a linear term B + Cr; (3) when volume is fixed, the parts of W and U that come from 
the linear term are almost constant. Our initial and simpler explanation of strong WU correlations (Ref. [I) was based on the 
dominance of close encounters, i.e., that it is only the nature of the repulsive part of the potential that matters for the strong 
correlations. This explanation, however, is adequate at high pressure / density only. It does not explain the requirement of fixed 
volume, nor the fact that strong correlation is observed even at zero pressure, as well as for the low-temperature / low-pressure 
(classical) crystal. To see that an explanation at the individual pair interaction level is generally inadequate, consider Fig. [TJb) 
which shows a scatter plot of single-particle energy and virial. These are sums over the pair interactions a given particle has 
with its neighbors; summing over all particles gives the total potential energy and virial, respectively. If strong correlation held 
at the level of single pair interactions, it would also hold at the particle level, but it clearly does not. This emphasizes that strong 
correlation is a collective effect, as detailed in Paper 11. 

In this paper we elaborate on the statistical mechanics and thermodynamics of strongly correlating liquids, and present re- 
sults from computer simulations showing that strong virial /potential energy correlations are present even in non-equilibrium 
processes. The purpose is to present a number of new results supplementing those of Paper II in order to broadly illuminate the 
properties of the class of strongly correlating liquids. Together Papers II and III give a fairly complete characterization of the 
properties of a strongly correlating liquid at one state point, as well as at different state points with same density. Paper IV in 
this serieJSgoes on to consider varying-density curves of "isomorphic" state points in the state diagram, which are characterized 
by several invariants; such curves exist only for strongly correlating liquids. 

The organization of this paper is as follows. Section|ll]begins with a discussion of the scaling properties of systems with inverse 
power-law (IPL) potentials, the natural starting point for a discussion of the hidden scaling properties of strongly correlating 
liquids. This is followed by a generalization to allow an extra term depending on volume only. Some, but not all, scaling 
properties of IPL systems are inherited by this generalization. Following this, in Sec. Ill we discuss the "extended inverse-power 
law"(eIPL) potential introduced in Paper II, which includes the above-mentioned linear term. We illustrate with simulation 
results the key property that the linear term contributes significantly to the virial and potential energy fluctuations when the 
volume may fluctuate, but little when it is fixed. Hence it gives rise, approximately, to a volume-dependent term in the free 
energy of the type discussed in the previous subsection. This leads to an inherited approximate scaling property, which we refer 
to as "hidden scale invariance" since it is not immediately obvious from the intermolecular potential. The argument about how 
and why hidden scale invariance causes strong WU correlations makes no assumption about equilibrium. To emphasize this 
point. Sec. [iVjpresents results from non-equilibrium computer simulations of strongly correlating molecular liquids, in particular 
aging following a temperature jump, and crystallization, both at constant volume. The property of strong correlation is shown to 
apply even in these out-of-equilibrium situations. Section |V]discusses ensemble-dependence, in particular it is here shown that 
the virial / potential energy correlation is always stronger in the NVT ensemble than in the NVE one. The last main section. Sec. 
|VI| comprises two topics under the heading "thermodynamics of strongly correlating liquids". First we discuss the relation of 
pressure-energy correlation to the thermodynamic Griineisen parameter jq, showing that the slope 7 (Eq. (j6])) is larger than jq 
by roughly a factor involving the ratio of excess (configurational) to total specific heats (constant volume). This ratio is around 
two for many simple liquids.!^ The second part formulates the property of strong correlation in the energy -bond language known 
as "network thermodynamics".G2l 
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FIG. 1: (a) Scatter plot of total virial and potential energy (in Lennard- Jones units) for the standard single-component LJ liquid at T = 80K 
(Argon units) and near-zero pressure, simulated at constant volume (density p = 34.6 mol/1, Argon units, left panel) and constant pressure 
(1.5 MPa, Argon units, right panel), (b) Scatter plot of single-particle virial and potential energy for the same simulation as in the left panel of 
(a). The single-particle WU correlation is much weaker, R = 0.63, showing that collective effects are crucial for the correlation. 



n. PROPERTIES OF INVERSE POWER-LAW SYSTEMS AND GENERALIZATIONS 

The purpose of this section is to summarize the properties of inverse-power law potentials and identify which of these proper- 
ties are inherited by strongly correlating liquids and which are not. 



A. Inverse power-law potentials 



Inverse power-law (IPL) potentials - sometimes referred to as soft-sphere potentials - have been used in liquid state theory for 
many years as convenient model systems." Such potentials have a number of simple properties. IPL potentials 

have, however, been considered unrealistic because their predicted equation of state is quite wrong and because they have no 
stable low-pressure liquid phase and no van der Waals loop, problems which derive from the fact that IPL potentials are purely 
repulsive. Moreover, the correct IPL exponent fitting the Lennard-Jones (LJ) liquid is around 18 (Papers I and II, Refs. .16.19.20i i. 
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not 12 as one might naively guess from the repulsive r~^^ term of the LJ potential; this may have confused people searching 
from an effective IPL description of the LJ liquid. A major point made in this series of papers is that when interpreted correctly, 
IPL potentials are much more realistic than generally thought, because they describe well a number of properties of strongly 
correlating liquids. For reference we now briefly summarize the since long well-established properties of IPL liquids. 

Consider N identical particles in volume V interacting by a pair potential of the form v{r) — Ar^"^; we make the pair 
assumption for simplicity but note that the below argument generalizes immediately to any potential that is an Euler homoge- 
neous function of the position coordinates. From the standard partition function the (Helmholtz) free energy F is conveniently 
written''^' as the ideal gas term plus the nontrivial "excess" free energy, F — F^^ + F^.^. The first term is the free energy of 
an ideal gas at same volume and temperature, Fjd = ~NkBTh\{pK^) where p = N /V i s the particle number density and 
A = h/ \/2iTmkBT is the thermal de Broglie wavelength. The excess free energy is giverP^by 



Whenever n > 3 this expression leads to a free energy with a well-defined extensive thermodynamic limit (N oo)! " l '^ l 

It follows from Eq. (|7]l that the excess free energy of an IPL liquid is given as follows in terms of a function of density p to 
the power n/3 over temperature T, 0(p"/'^/r) (Klein's theoremii-^: 



Fex,iPL = NkeTcj) ■ (8) 

This implies that a number of derived quantities are also functions of p^^'^/T. As important examples, recall the following 
standard identities: The excess entropy: 5cx = —{dFcx/OT)v, the potential energy: U = Fcx + TSq^, the virial W — 
—V{dFcx/dV)T, the excess isothermal bulk modulus: = V{d^Fcx/dV^)T, the excess isochoric specific heat per unit 
volume: = -(T/y)(92i^ex/9T2)y, the excess pressure coefficient: f3^ = {l/V){dW/dT)v = -d^F^JdTdV. From 
Eq. ([sjl it follows that these three quantities are functions of the single variable p^^'^ /T; more accurately one has [where 
h{x) = x^'{x) - f2{x) - x(t>'{x), h{x) - {n/ifx^"{x) + [(vi/S) + {n/?,f]x^'{x), and U{x) = -^^^"(a;)] 



5ex,iPL = Nks h [p'^'Vt) , (9) 
C/iPL - NkBTf2 [p"^yT^ , (10) 



W^iPL = '^NkBTh[p''^yT'j , (11) 

^T>L = pfcBr/3(p"/Vr) , (12) 
4Vl = pksh , (13) 

PvW = '^pkB h {p'''Vt) . (14) 

The functions /i, /4 all depend on n, but for simplicity of notation we have not indicated this explicitly. Dividing across 
by the dimensional factors on the right hand side (for example fc^T in the case of potential energy and virial), one arrives at 
dimensionless forms of the excess entropy, potential energy, etc. that are functions of p'^/^ /T only. 

Turning now to the dynamics, consider the standard molecular dynamics (MD) case where the equations of motion are 
Newton's equations. Suppose ri(i) (i = 1, A^) is a solution to Newton's equations. Then it is straightforward to show that 

v'^p {t) — aVi{\t) is also a solution if a^^"'^'^^ = A^. In particular, if Vi{t) refers to equilibrium (NVE or NVT) dynamics at a 

state point with density po and temperature Tg, then v'f'\t) = ar,;(Ai) refers to equilibrium dynamics at density pi — pa/a^ 
at temperature Ti = Toa^X^ (temperature scales as the mean-square velocity and velocities get a factor aX). Using the above 
relation between a and A this implies 



Ti = a-"To = ^ To . (15) 

This means that two states with different densities and temperatures but same p^^^^/T have dynamics that scale into one 
another by simple scalings of space and time. In particular, if for any quantity A one defines the relaxation time ta via 
TA = f^{A{0)A{t))dt/ (A^), it follows that any two states with same p'"'^^/T have same "reduced" (dimensionless) relaxation 
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time ta, if this quantity is defined by ta = TA/to where the characteristic thermal time to is defined by to = P ^^"^ \frnjkj^ . 
Similarly, if one defines the reduced diffusion constant Dhy D = D/Dq where -Do = p~^^^/to = p~^l^ ^fkj^Tjm, then t) is 
the same for the two states. Summarizing, 

T = h (f'^lT) , (16) 

D = h {p"^'/t) . (17) 
Finally we note that it foUows from the above scaling property that 

^jPL(rl^\...,r^^) ^ ?7/PL(ri,...,rjv) ^^^^ 

Thus the Boltzmann factors of the two configurations are the same and, consequently, the scaling of the dynamics holds also 
for stochastic dynamics; this observation, in a generalized form, is the starting point of Paper IV in this series. By the same 
argument the structure of states with same p"/^/T are identical, provided lengths are scaled by p~^l^. 



B. Inheritance of scaling properties by generalized IPL potentials 

This section discusses how a large class of potentials inherits a number of IPL properties to a good approximation, thus 
justifying the term "hidden scale invariance". Consider a general potential between particles i and j, rewriting it (as can always 
be done) as a sum of an IPL potential plus the difference, denoted "diff": 

+<f^i3)- (19) 

For any configuration, (ri, rjv), the potential energy is then the sum of an "IPL" term and a "diff" term, and the excess free 
energy is given by 

g-Fex/fcsT ^ y" ^ ^^g-(7iPL(ri,...,rjv)/feBTg-i7ditt(ri,...,rjv)//sBT _ ^-20) 

We now investigate consequences of the assumption (Sec. Ill) that C/dift to a good approximation is only a function of volimie: 
{7(jiff(ri, r^v) = /(V^) - at least for states that carry Boltzmann weights of any significance. The approximate identity 
C^diff (ri, rjv) = /(^) means that the second exponential can be moved outside the integral, and we get: 

g-Fex/fesT ^ g-/(V)/feBT /■ ^^_ J^g-(7iPL(ri,...,r^,)/feBT_ fr^y^ 



From this follows directly that 



^ ex 



J{y) + -Fex.iPL = !iy) + NkBT<i> (fl^lT) , (22) 



which implies 



5ex = 5ex,IPL = A^fcB/l , (23) 

V = /(F) + C/iPL = /(V) + NksTh {/''^IT) , (24) 

W = -r{y)V + = -r{V)V + \NkBTh {p^'Vt) , (25) 

i^- = Vf"{V) + Jfl^ipL = Vf"{V) + pksTh if'^T) , (26) 

= cf?iPL = pkBh (p^^Vt) , ill) 

P^v = Pv,iPL = IpkB h [p'^'^It] . (28) 
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While the systems under consideration here have the same excess entropy as the "hidden" IPL, several quantities have con- 
tributions from the /(F)-term. These quantities do not obey IPL scaling. In contrast, the scaling behavior for dynamics and 
structure is inherited: Consider two state points with the same p"/^/T. For the pure IPL system (f{V) — 0) the two state points 
have the same dynamics and structure as argued in the previous section. Letting f{V) ^ simply shifts the energy surface, 
which changes neither the dynamics nor the structure. This scenario - scaling of the dynamics, but not the equation of state - 
is exactly what is observed experimentally for a large number of viscous liquid s. For example in van der Waals liquids relaxation 
times are found to be a function of p"-/^/T (using n as an empirical parameter) the scaling does not apply to the (excess) 

pressure with the exponent determined from the scaling of relaxation time, as required for IPL scaling.^"* 

In the following section we provide numerical evidence that there are indeed systems that to a good approximation fulfil the 
assumption introduced above that Uditf is a function of volume only. 



III. LENNARD- JONES AS A GENERALIZED IPL POTENTIAL: THE EXTENDED INVERSE POWER-LAW POTENTIAL 

APPROXIMATION 



In this section we examine the extent to which the Lennard- Jones (LJ) potential may be approximated by an "extended inverse 
power-law" potential (including a linear term, Eq. ( [33] l below), by considering the fluctuations at a particular state point of the 
LJ fluid. We choose a state point whose pressure is near zero, because here it is particularly clear that single-pair effects are 
insufficient to explain the strong WU correlation (Fig. 1). 

The analysis of Paper II took its starting point in assuming that the approximating inverse power law should match the potential 
closely at a particular value of interparticle separation. An important conclusion of the analysis was, however, that the success 
of the IPL approximation derives not from its behavior near any particular r-value, but rather from the fact the difference from 
the real potential is close to linear over the whole first-peak region. In this section we re-examine the idea that the fluctuations 
are well described by an inverse power-law potential and the argument for why the difference term almost does not fluctuate. 
We show explicitly that the latter contributes little to the fluctuations at constant volume, but significantly when the volume is 
allowed to fluctuate as in the NpT ensemble. This demonstrates that the LJ potential is of the type considered in the previous 
subsection. 

We wish to determine to what extent the LJ potential 

VLjir) = 4e ((a/r)i2 - (a/r)^) (29) 

can be matched, for the purpose of describing fluctuations of potential energy and virial at fixed volume, by an inverse power 
law 



w/Pi(r) = Ar-"[+B]. (30) 

Here B indicates an optional constant. To start with, how should the exponent n and the coefficients A and B be chosen? An 
obvious choice, followed in the first part of Paper II, is to require that the two potentials, vlj and vjpl, should agree as much as 
possible around a particular value of r, denoted tq. Given ro, if we require the functions and their first two derivatives to match 
at rg, this determines all three parameters A, n and B. The exponent n is given (Paper II) by 



n = nW(.o)^-^^^^^-l. (31) 

Here the notation n^^\ro) refers to one kind of r-dependent effective inverse power-law exponent, based on the ratio of the 
second and first derivatives. For vipL{r) this simply returns n. Otherwise it gives a local matching of the vjpL{r) to VLj{r). 
This leaves effectively one parameter to vary, namely tq, which must be less than the minimum r™ — l^l^a where n^^' diverges. 
The parameter tq may be chosen to optimize the match of the fluctuations in total energy and virial. For an NVT simulation at 
T — 80K and near-zero pressure, the best fit was obtained with n = 19.2 (while the exponent implied by the slope [Eq. (j6]l], 
7 — 6.3, was slightly smaller, 18.9). 

Later in Paper II it was demonstrated that there is no particular reason why the potentials should match close at a particular 
value of r, since the fluctuations have contributions from the whole first-peak region, including beyond the potential minimum. 
The reason that any kind of matching is possible over this region - where ulj clearly does not resemble a power law - is that 
a linear term may be added to the power-law potential almost without affecting the fluctuations as long as the volume is held 
constant. The analysis of Paper II, which also included an in-depth treatment of the perfect LJ (fee) crystal which is also strongly 
correlating, showed that the more relevant r-dependent effective exponent is the higher order n*^^^ defined by 
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FIG. 2: (Color) Comparison of g{r) for simulations using the Lennard- Jones potential and two inverse power-law potentials: the r^^^ 
repulsive term in VL.jir) and the inverse power-law potential that optimizes the agreement in the fluctuations of potential energy and virial by 
minimizing Eq. \i6) . The left panel shows these at density 0.82 and temperature 0.67 (dimensionless units), the right one at density 0.90 and 
temperature 0.80 (where the r~^^ potential leads to crystallization). 



This also returns n for vjpl {r), but since it does not involve the first derivative, it returns n even if a linear function of r is added 
to the potential as in the extended inverse power-law potential (eIPL) defined by 



VeiPL = Ar-"" + B + Cr . (33) 

This potential fits the LJ potential very well around its minimum (Paper II) and thus includes part of its attractive part. 

There are several possible ways of choosing the "best" eIPL to match the real potential. These will give slightly different 
exponents and coefficients A, B and C. We do not investigate them here; rather the purpose is to validate the basic idea of 
the extended inverse power-law (eIPL) approximation. Therefore we choose a simple matching scheme, whereby we match 
the fluctuations to those of the inverse power-law potential, without including a linear term, in order to determine n and A. 
For simplicity we take the exponent directly from the observed fluctuations: n = 87 where 7 is defined in Eq. (j6]). To fix the 
coefficient A, agreement with the potential energy and virial fluctuations is optimized by proceeding as follows: For a given 
configuration generated in an LJ molecular dynamics simulation we calculate the LJ potential energy Ulj and the power-law 
potential energy Ujpl, similarly the corresponding virials W^j and Wjpl. The difference quantities C/dift and Wdiff are defined 
as 



(7diff = Ulj - Uipl , (34) 
Wdiff = Wlj - WiPL ■ (35) 

A perfect match of the fluctuations would mean that C/diff and Wdiff have zero variance. Therefore we choose A to minimize the 
sum of the relative "diff" variances: 



((A£/d.ff)^) ((Al^diff)^) 

For the near-zero pressure state point used in Fig. 1 of Paper I the exponent determined from 7 is n = 87 = 18.9 and the 
optimal value of A is 1.3437e(T". Before examining the difference potential, what do we get if we simulate with the matched 
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FIG. 3: (Color) (a) Illustration of the difference between the Lennard-Iones (LJ) potential VLjir), the empirically matched inverse power-law 
potential vipLij-) with A = 1.3437e(7" and n = 18.9, and their difference VMsij-). (b) Linear fit, «iin(r) = min(0, —3.6635 + 2.2756r/cr), 
to Wdiff between 0.95(t and 1.5cr, and the remainder i'rcst(r) = Udiff(f') — v\m{r) (full black curve). 



inverse power-law potential? Figure |2] shows the radial distributions g{r) obtained for the above state point and another with a 
higher density and temperature, for three potentials: LJ, the repulsive r^^^ term of the LJ potential, and the optimal IPL potential 
with n = 18.9. We used the same inverse power-law potential at both state points (i.e., we did not adjust A and n to match the 
second state point). The first thing to note is that the n = 18.9 potential gives a structure much closer to that of LJ than does the 
repulsive n = 12 term alone, in particular the latter system has crystallized at the higher density and temperature. The second 
point is that there is still a difference between the LJ and the n = 18.9 IPL, present in both state points. The first peak in the LJ 
system is slightly higher and narrower, although its position is barely altered. Thus the real potential gives a slight increase of 
order - the difference in coordination number is less than 0.1 (integrating to the first minimum after the peak). 

Figure [3| a) shows the LJ potential, the IPL potential with parameters optimized as described above, their difference, and the 
radial distribution function. As shown in Fig.[3|b), the main part of the difference potential v^is{r) is nearly linear Thus a good 
approximation to the real potential is the eIPL potential of Eq. p3] l for r less than a cut-off rc, and zero otherwise. Neglecting 
the small value of vjpL{rc) ^ lO^^e, the cut-off is given by —B/C. For the fit shown in Fig.[3](b), ~ l.Glcr. 

What are the implications of the linear term? A linear term in the pair potential contributes a term proportional to the sum 
of all bond lengths to the total potential energy. It was shown in Paper II that at constant volume this sum is a constant in one 
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FIG. 4: (Color) Effect on fixed-volume fluctuations of adding a linear term to the inverse power-law (IPL) potential. The linear term is that 
shown in Fig. [3|b). Configurations were generated by an NVT simulation using the LJ potential, and the different determinations of energy 
(LJ, IPL and eIPL) and virial were computed on these configurations. The dashed lines indicate a perfect match. Including the linear term 
when computing the energy improves the match to the true (LJ) fluctuations (the correlation coefficient goes from 0.950 to 0.970), while it 
reduces the match to the virial (the correlation coefficient goes from 0.987 to 0.971, which is probably related to the fact that the pair virial is 
discontinuous at Vc - we thus find that smoothing the linear part around Vc restores the match somewhat). The insets show the pair potentials 
and virials: Brown dashed lines: LJ, black lines: IPL, red lines: elPL. - The overall conclusion from Fig. |4]is that the addition of the linear 
term induces little change in the fluctuations. 



TABLE I: Variances of potential energy U and virial W , and of various contributions to U and W, of two different ensembles at the LJ state 
point given by p = 0.82 and T — 0.67 (dimensionless units). 



Ensemble 


Quantity 


LJ 


IPL 


diff 


lin 


rest 


NVT 


U 


0.0231 


0.0225 


0.0075 


0.0085 


0.0063 




W 


0.1468 


0.1402 


0.0227 


0.0301 


0.0350 


NpT 


u 


0.0484 


0.0320 


0.0665 


0.0539 


0.0144 




w 


0.1704 


0.1997 


0.0589 


0.0417 


0.0430 



dimension, and it was argued that it is approximately constant in three dimensions. The difference is because of two things: 
First, in three dimensions there are contributions to bond-length changes from transverse components of relative displacements 
between the two particles defining a bond; secondly, within the eIPL the potential is only linear up to Vc and, moreover, significant 
deviations of w^iff from linearity occur at r's smaller than r^)- 

The argument that the linear term contributes little to the fluctuations depends on all bond lengths being less than rc- In one 
dimension, at moderate temperatures, a single-component system has a rather well-defined nearest-neighbor distance, which at 
densities where the pressure is not too negative will be less than Vc- In a three-dimensional liquid, however, the nearest-neighbor 
distance is not as well-defined - the radial distribution function does not go to zero after the first peak. Therefore there will 
always be fluctuations at Vc as the lengths of bonds fluctuate back and forth across rc, so the sum of bond lengths which are less 
than Tc will fluctuate. In Paper II it was shown that for a three-dimensional (classical) crystal at low temperature - where this 
is not an issue because g{r) does go to zero after the peak - the correlation coefficient R becomes very high, over 99.5%, as 
(but not 100%). 

We can check directly the effect of adding the linear term to the inverse power law. Figure]?] shows scatter plots of IPL (black) 
and eIPL (red) energy and virial plotted against the true LJ values. These were calculated for a set of configurations drawn from 
an NVT simulation using the true (LJ) potential. Including the linear term makes little difference. It somewhat improves the 
match to the energies, though not to the virials (possibly due to the discontinuity in the pair-virial at rc). The WU correlation 
coefficient of the eIPL potential energy and virial is 0.917 (compare to the true (LJ) value of 0.938 and the pure IPL value of 
1.0). 

It is instiTictive to repeat the above for configurations drawn from an NpT simulation at the same state point, i.e., with pressure 
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FIG. 5: (Color) Comparison of potential energy calculated using VLjir), vipL(r) (black points) and VeiPLij-) (red points) for configurations 
drawn from an NpT simulation using the LJ potential. The potential energy, in particular, is very poorly represented by the power-law 
contribution when the volume is allowed to fluctuate. In fact the correlation between Uipl and Ul.j is not only weak, it is negative. Including 
the linear term makes a huge difference here, yielding a correlation coefficient of 0.977 between UeiPL and Ulj (changed from -0.201). The 
slope is somewhat less than unity, indicating that there are significant contributions from pair distances beyond Tc (i.e., from the "rest" part of 
the potential). The linear term affects the virial fluctuations much less, presumably because the derivative of the potential is dominated by the 
IPL term. 



chosen as the average pressure in the NVT case. The results are shown in Fig. [5] Here it is clear that the IPL potential represents 
the potential energy fluctuations very poorly (black points), while adding the linear term makes a substantial difference (red 
points). As in the NVT case the linear term affects the virial fluctuations much less. This is presumably because when taking 
the derivative to form the virial, the IPL term gets multiplied by n = 18 while the linear term gets multiplied by minus one, and 
is thus reduced considerably in significance (compare the insets of Fig. [4]). 

The size of the variances of the different terms are compared in Table |l] We do not make a detailed analysis of the variance 
(taking into account cross-correlations, etc). In the NVT case the IPL contributions are of similar size to the full (LJ) fluctuations 
- naturally since we explicitly optimized this - and the "diff" contributions are small compared to the IPL ones. In the NpT 
case, on the other hand, the "diff" contributions to the fluctuations of U are more than double the IPL ones; this is not the case 
for the "diff" contributions to W, though they are still a larger fraction of the total than in the NVT case. These numbers are 
consistent with Fig. [5] The "diff" contributions to the energy must be larger than the IPL ones because the latter are negatively 
correlated with the true energy. The fact that the variance of J/diff is smaller than the sum of those of U\m and J/rest, in the NVT 
case, indicates that the latter are negatively correlated. This is presumably due to bond-lengths around which alternately are 
counted as part of U\in and as part of t/rest- This effect is less noticeable in the NpT case; there we see clearly that fluctuations in 
U\ia account for most of those in C/,est- 

Based on the above we can now answer the question; Is it possible to predict whether or not a liquid is strongly correlating by 
inspection of its potential (i.e., without simulating virial and potential energy fluctuations)? For liquids with particles interacting 
by pair potentials, the answer is the affirmative: The liquid is strongly correlating if the potential around the first peak of the 
structure factor (the typical interparticle distance) may be fitted well by an extended inverse power law. For more general 
potentials, the situation is more complex. Thus it is possible to construct many-body potentials with angular dependencies 
which scale like inverse power-law pair potentials; these have 100% WU correlation because this property follows whenever 
the potential is an Euler homogeneous function. In most realistic cases, however, systems with angular dependencies are not 
expected to be strongly correlating; likewise potentials with two length scales will generally not be strongly correlating. An 
example of the former is the coarse-grained model of water using a short-range many-body potential recently introduced by 
Molinero and Moord^ that reproduces water's properties with surprising accuracy. This potential is not strongly correlating, 
because close to water's density maximum the virial / potential energy correlation coefficient R must be close to zero (Paper 1). 
Examples of potentials with two length scales, that for this reason are not strongly correlating, are the Jagla potential^ and the 
Dzugutov potential.'^ Likewise, the addition of Coulomb terms to an LJ-type liquid generally ruins strong correlations (Paper I). 
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FIG. 6: (Color) Computer simulations of virial and potential energy during the aging of two strongly correlating molecular liquids following 
temperature down-jumps at constant volume (NVT simulations), (a) The asymmetric dumbbell model at density p — 1.109 g/ml. The liquid 
was first equilibrated at T=300 K. Here simultaneous values of virial and potential energy are plotted at several times producing the green 
ellipse, the elongation of which directly reflects the strong WU correlation in equilibrium. Temperature was then changed to T=240 K where 
the red ellipse marks the equilibrium fluctuations. The aging process itself is given by the blue points. These points follow the line defined by 
the two equilibrium simulations, showing that virial and potential energy correlate also out of equilibrium, (b) Similar temperature down-jump 
simulation of the Lewis-Wahnstrom OTP systemP^The colors have the same meaning as in (a): Green marks the high-temperature equilibrium 
(T=600 K), red the low-temperature equilibrium (T=450 K), and blue the aging towards equilibrium. - In both (a) and (b) the slope of the 
dashed line is not precisely the number 7 of Eq. l|6l because the liquids are not perfectly correlating; the line slope is (At/ AVI^) /{(AC/)^) , 
see Paper I, a number that is close to 7 whenever the liquid is strongly correlating, (c) Virial and potential energy for the asymmetric dumbbell 
model as functions of time after the temperature jump of (a); in the lower subfigure data were averaged over 10 ps. Virial and potential energy 
clearly correlate closely. 



IV. OUT-OF-EQUILIBRIUM DYNAMICS IN MOLECULAR MODELS 

According to the extended inverse power-law (eIPL) explanation detailed in the previous section, strong WU correlations 
characterize all configurations of the LJ liquid at a given volume. This means that the correlations should be there also under 
non-equilibrium conditions if the volume is kept constant. In this section we present numerical evidence that this prediction is 
indeed fulfilled, even for molecular models, provided they are strongly correlating in their equilibrium WU correlations. 



A. Temperature down-jump simulations of three molecular model liquids 



Figure [6|a) shows the results for a temperature down-jump at constant volume, starting and ending in equilibrium (NVT 
simulations).^'' The system studied is an asymmetric dumbbell liquid consisting of two different-sized LJ particles glued together 
by a bond of fixed length, with parameters chosen to mimic toluene.l^The system was first equilibrated at 300 K. The green ellipse 
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FIG. 7: (Color) Virial versus potential energy after a temperature down-jump at constant volume applied to SPC water, which is not strongly 
correlating (colors as in Fig. |6](. (a) SPC water at 1 atm equilibrated at T=260K, subsequently subjected to an isochoric temperature down 
jump to T=160 K. Clearly W and U are not strongly correlated during the aging process, (b) Same procedure starting from a 2GPa state point. 



consists of several simultaneous instantaneous values of U and W in equilibrium at T=300 K. The strong W\J correlation is 
revealed by the elongation of the ellipse (i? = 0.97; 7 — 6.1). When the liquid is similarly equilibrated at 240 K, the red blob 
appears. To test for correlation in an out-of-equilibrium situation we changed temperature abruptly from the 300 K equilibrium 
situation to 240 K. The blue points show how virial and potential energy evolve following the temperature down jump. Clearly, 
strong Wl] correlations are present also during the aging towards equilibrium. Figure|6|b) shows the same phenomenon for the 
Lewis-Wahnstrom ortho-terphenyl (LW OTP) model which consists of three LJ spheres at fixed length and angle with parameters 
optimized to mimic ortho-terphenyl.^'' This liquid is also strongly correlating (i? = 0.91; 7 = 7.6). The colors are as in Fig. 
[6|a): Green gives a high-temperature equilibrium state (T = 600 K), red a low-temperature equilibrium state (T = 450 K), 
and the blue points show the aging towards equilibrium after changing temperature from 600 K to 450 K. The picture is the 
same as in Fig. |6|a): The blue points follow the dashed line. Thus virial and potential energy correlate strongly also for far- 
from-equilibrium states. Figure |6](c) plots W(i) and U (t) after the temperature jump for the Fig. |6|a) data for the asymmetric 
dumbbell liquid. W{t) and U (tjfollow each other closely on the picosecond time scale as well as in their slow, overall drift to 
equilibrium. 

What happens when the same simulation scheme is applied to a liquid that is not strongly correlating? An example is SPC 
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FIG. 8: (Color) Crystallization of the supercooled Lewis-Wahnstrom ortho-terphenyl (OTP) liquid where each molecule consists of three 
Lennard-Jones spheres with fixed bond lengths and angles.^' (a) Pressure (right) and energy (left) monitored as functions of time during 
crystallization at constant volume. Both quantities were averaged over 1 ns; on this time scale the pressure / energy fluctuations directly reflect 
the virial / potential energy fluctuations. The horizontal dashed lines indicate the liquid (upper line) and crystal (lower line), the averages 
of which were obtained from the simulation by averaging over times 0-2 /is and 5-10 fis, respectively. Both liquid and crystal show strong 
correlations, and the correlations are also present during crystallization. Inset: Crystal structure from the simulation, (b) Radial distribution 
functions of liquid and crystalline phases. The two spikes present in both phases come from the fixed bond lengths. 

water, where the hydrogen bonds are mimicked by Coulomb interactions.'^ Figure [t] shows results of simulations of SPC water 
at two different densities, (a) corresponding to low pressure and (b) to very high pressure. In the first case virial and potential 
energy are virtually uncorrelated {R = 0.05, T = 260 K); in the second case correlations are somewhat stronger {R — 0.34, 
T = 260 K), though still weak. As in Fig. |6] green denotes the initial, high-temperature equilibrium, red the low-temperature 
equilibrium, and blue the aging towards equilibrium. Clearly, for this system W and U are not closely linked to one another 
during the relaxation towards equilibrium. 

B. Pressure and energy monitored during crystallization of a supercooled liquid: The Lewis-Wahnstrom OTP model 

A different far-out-of-equilibrium situation is that of crystallization of a supercooled liquid monitored at fixed volume and 
temperature. To the best of our knowledge crystallization of the LW OTP model has not been reported before, but Fig. |8]shows 
that for simulations over microseconds the supercooled liquid crystallizes at T=375 K and p = 1.135g/cm'^. In the crystal the 
LJ spheres are aiTanged in a slightly deformed bcc lattice where the molecules have otherwise random orientation. The crystal 
is shown in the inset of Fig. |8ja). Figure |8ja) shows how time-averaged pressure and energy develop during crystallization. 
Contributions to pressure and energy from the momentum degrees of freedom are virtually constant after averaging over 1 ns, 
so strong WU correlations manifest themselves in strong averaged pressure /averaged energy correlations. Clearly averaged 
pressure and averaged energy follow each other closely also during crystallization. This confirms the above finding, as well as 
those of Ref. that strong correlations apply also for the crystalline phase of a strongly coiTelating liquid. Note that the slope 7 
is virtually unaffected by the crystallization. The persistence of strong virial / potential energy coiTelation during crystallization 
- and the insignificant change of 7 - are noticeable, because physical characteristics are rarely unaffected by a first-order phase 
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FIG. 9: (Color) WU plot for the asymmetric dumbbell model for various states at the same volume. The upper right comer shows data for 
simultaneous values of virial and potential energy for four equilibrium simulations (T=240-350 K). When quenching each of these to zero 
temperature in order to identify the inherent states, the crosses are arrived at. The intermediate points are glasses prepared by different cooling 
rates: out-of-equilibrium systems generated by cooling in 1 ns from 240 K to the temperature in question. This plot shows that strong virial / 
potential energy correlations are not limited to thermal equilibrium situations. 

transition. These simulations show that the property of strong virial / potential energy correlations pertains to the intermolecular 
potential, not to the particular configurations under study. Figure [Hj^b) shows the radial distribution functions for the liquid and 
crystalline phases. 

C. Glasses and inherent states 

The above out-of-equilibrium simulations show that the property of strong WU correlation is not confined to thermal equilib- 
rium. This is consistent with the elPL description given in the previous section, but was shown here to apply even for molecular 
models with fixed bonds. It appears that strongly correlating liquids have a particularly simple configuration space. These re- 
sults have significance for any out-of-equilibrium situation. Consider the potential energy landscape picture of viscous liquid 
dynamics^' ^- according to which each configuration has an underlying inherent state defined via a deepest-descent quench, 
a stat e that con tains most information relevant to the slow dynamics, which may be regarded as jumps between different inherent 
stateJ^fSlillll Figure[9]shows a WU plot of the asymmetric dumbbell model in different situations at same density: equilibrium 
states (upper right) and their corresponding inherent states (lower left, one quench per temperature), and glasses at different tem- 
peratures in between. A glass is an out-of-equilibrium state, of course, and inherent states may be regarded as zero-temperature 
glasses. The plot shows, once again, that strong correlations are present also far from equilibrium. 

V. ENSEMBLE DEPENDENCE OF THE CORRELATION COEFFICIENT 

Most of the simulations of Paper 1 were done in the NVT ensemble. An obvious question is how the correlations differ 
between ensembles. This section develops the necessary theory needed to answer this question and compares first the NVT and 
NVE ensemble, then the NVT and NpT ensembles. 

It is well known that although simple thermodynamic averages are independent of ensemble,^^ fluctuations are generally 
ensemble dependent. Consider two ensembles, one with extensive variable F held fixed, and one with its conjugate intensive 
variable / held fixed (defined so their product is dimensionless). The other parameters defining the ensembles are the same. The 
covariance of observables A and B in the two ensembles are related aJ^^^H 
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{AAAB)j. = {AAAB)f + [ ) . (37) 



OF J \df' 'V \9f 
A. NVT versus NVE 

To compare the NVT and NVE ensembles we take F and / as the energy E and the inverse temperature /? = \/{kBT), 
respectively, keeping the volume fixed in both cases. Noting that d/dp — —ksT^d/dT and {d[3/dE)v — l/(— fcsT^Ci/) 
where Cv is the extensive isochoric specific heat, Cy = Vcy, the covariance of U and W is given by 

/A7-rAw\ /A7-rAw\ ksT"^ d{U)NVT d{W)NVT 

{AUAW)nve = {AUAW)nvt - —7^ (38) 

= {AUAW)nvt - kBT^C'^V(3'^/Cv • (39) 

Here we have introduced the excess parts of the isochoric specific heat and pressure coefficient, Cf? and (3y, respectively (note 
the change in notation from Paper 11 where the superscript "conf" was used - "ex" is however more standard in liquid state 
theory). These two quantities are given by 

C%- = I ^ I = Cv - ^Nks (40) 
I^^A--^^ (4.) 
For the variances one has 

{iAU)^)NVE = {{AUf)NVT - kBT\C-^f/Cv (42) 
{{AWf)NVE = {{AWf)NVT - kBT\Vl3'^f/Cv ■ (43) 

The above implies that the NVE W [/-correlation coefficient Re (in the following subscripts E or T indicate the NVE or NVT 
ensemble, respectively) is given by 

^ {AUAW)NVT-kBT^C'^VP'^/Cv 

^{{au)^)mvt - kBT^c-^)ycvy'{{Awy)NVT - kBT^ivp-^Y/Cv ■ 

We wish to express the right side in terms of the NVT coefficient Rt- The definition of Rt implies 





(d{U)\ 


V — 


V dT )y 


ex 


(d{W/V) 


V — 


\ dT 



, {AUAW)%y^ 

HAW) )wT - ^2 ^(^f,).^^^^ ■ (45) 



Inserting this into Eq. (44 1 and making use of the fluctuation relations {AUAW)nvt — kBT'^VPy and {{AUY) nvt 



kBT^Cy (see, e.g., appendix B of Paper 1) gives 



^ ksT^Vp'^ - kBT^C'^VP'^/Cy ^^^^ 

After squaring and cancelling factors of fc^T^ and V(3y we get an expression relating Re to Rt, 

- II Rl C^^/Cv ■ ^ ^ 
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system 

LJ 

LJ 

LJ 

LJ 
KABLJ 



1.00 
1.00 
0.82 
0.82 
1.2 



1.00 
0.80 
0.80 
0.67 
0.47 



0.991 
0.991 
0.943 
0.949 
0.936 



C^/N_ 
1.5 
1.7 
0.90 
1.3 
2.1 



Eg. ([47 I 
0.982 
0.981 
0.912 
0.909 
0.859 



Re 
0.983 
0.981 
0.918 
0.904 
0.862 



TABLE II: Check of relation lj47j between Rt and Re for the LJ and the Kob- Andersen binary Lennard- Jones (KABLJ) fluids. The units for 
p and T are the dimensionless units defined in terms of the length and energy parameters a and e for the interactions of the large particles. The 
excess isochoric heat capacity was calculated from the potential energy fluctuations in NVT ensemble. 



To get a feel for the relation, divide by R^. This yields 



Re 
Rt 



1 - C-^/Cy 
1 — R^Cy I Cy 



(48) 



First one notes that when Rt = 1, the denominator on the right side becomes equal to the numerator and Re = 1- That is, the 
property of perfect correlation is independent of (fixed-volume) ensemble. When Rt < 1, the denominator becomes greater 
than the numerator, and so R^ < R^. That is, the correlation coefficient is smaller in the NVE ensemble than in the NVT 
one. How can we understand this? Consider the set of WU points sampled by the system during an NVE trajectory; this is an 
elongated blob in the WU diagram. Changing the energy will cause the blob to move along a line almost parallel with the long 
axis of the blob."* Switching to the NVT ensemble is equivalent to superposing several of these collinear blobs on top of each 
other — the result is necessarily longer, but not wider This corresponds to higher correlation. 
Simulation data confirming relation (47 1 are presented in Table [ll] 



B. NVT versus NpT 



In the NpT ensemble, where volume is allowed to fluctuate, we must consider different variables. The natural variables to 



coiTelate are the excess enthalpy Hex = U + pV and the volume V. We use again Eq. ( 37 1, but now take FusV and / as p/3, 
the pressure times inverse temperature, keeping temperature constant. Eq. ( [37| becomes 



{AAAB)nvt = {AAAB) 



(49) 



The details of the calculation, which are somewhat tedious, are given in Appendix A. The result for the H^^V correlation 
coefficient is rather simple, though: 



R 



1 



H,^V,NpT 



(50) 



where 



and 



a = {AUAW) NVT + N{kBTf 



(51) 



62 ^ KTVkBT{{AUf 



INVT- 



(52) 



Notice that Rh^^v,NpT is strictly less than unity - even for perfectly correlating liquids (that is, with perfect WU correlations 
in the NVT and NVE ensembles). For the Lennard- Jones simulation of Fig. [T|(b) Rh^^v,NpT — 0.86 (recall the NVT W, U 
coiTelation coefficient for the same state point is 0.94). Unlike the situation when comparing the NVT and NVE ensembles, 
there does not seem to be a simple relation between the two correlation coefficients. It seems likely, though, that the HcyV 
correlation in the NpT ensemble is generally smaller than the WU correlation in the NVT ensemble; thus the property of strong 
correlation is less evident in the NpT ensemble. 
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VI. THERMODYNAMICS OF STRONGLY CORRELATING LIQUIDS 



The property of strong virial / potential energy correlation not just refers to microscopic properties that are only accessible in 
simulation, it also has consequences for the liquid's thermodynamics as well. The first subsection below relates the slope (Eq. 
(j6]l) to the Griineisen parameter, the second subsection shows how to give a general thermodynamic formulation of the property 
of strong correlations. 



A. Relation to the Griineisen parameter 

The Griineisen parameter was originally introduced to characterize the volume dependence of normal modes of a crystalP^I^Il 



where oji is the frequency of the ith normal mode. By assuming that ji is the same for all modes and denoting the common value 
by 7(5, one can derive the Mie-Griineisen equation of state,'^ 



P+di=^^^- ^^^^ 

Here p is pressure, u{v) with v — V/N is the "static" energy of the crystal per atom (the energy of the force-free configuration 
about which vibrational motion occurs), and i?vib is the vibrational energy. In general jq depends on volume, but this depen- 
dence is typically small enough that it can be neglected. From Eq. (|54]) it follows that, if E is the total, thermally averaged 
internal energy, one has {dp/dT)v = {-fG/V){dE^ii,/dT)v = {-fG/V){dE/dT)v, i.e.. 



This expression is the slope of the pressure versus energy curve at fixed volume, analogous to the 7 of Eq. (|6]l but for the presence 
of the kinetic terms (recall that for an IPL liquid 7 = n/3 — {dW/dU)v)- If ctp is the coefficient of thermal expansion, Kt 
the isothermal bulk modulus, and cy = Cy/V the isochoric specific heat per unit volume, Eq. ( [55] l implies via standard 
thermodynamic identities 



7g = — ■ (56) 

Cv 



This relation allows 7g to be determined from experimentally accessible quantities; in fact Eq. ( [56] l can be taken as a thermody- 
namic definition of 

There have been suggestions of how to connect the so-called density scaling exponenP^J25|_ the one controlling the relaxation 
time via the variable p'' /T - with the Griineisen parameter, notably by Roland and c o workers .l^^llolll] jn Ref.fW equality of jg 
and 7 was argued theoretically. More recently, Roland and Casalini showed"*' that equality is not consistent with experimental 
results; rather 7g is smaller than 7 by a factor of order three. This discrepancy was reconciled in the context of the entropy model 
for relaxation (by which the relaxation time is a unique function of the so-called configurational entropy 5c), by introducing a 
corrected Griineisen parameter defined via 



Here Acy is the difference of the isochoric specific heats per unit volume between the liquid and the glass. Because Acy 
is smaller than cy, one has 7™" > 70- By arguing from experimental data that the non-configurational part of the entropy 
(associated with vibrations, equal to the entropy of the glass) is independent of volume and assuming that ACy is constant, they 
derive 



Sc = ACy ln(TV^'^"') + const . 



(58) 
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Hence 7™" is the density-scaling exponent. 

We take a different approach to connecting the Griineisen parameter with the slope 7 (which provides a good estimate of the 
density scaling exponent, see Ref. 22 ). Instead of splitting the entropy, we split the pressure, into potential and kinetic parts and 
get from 7g = V{dp/dE)v 



^G = V^^^^-^^^ . (59) 



Expressing the temperature derivatives in terms of fluctuations (Appendix B of Paper I) gives 



{AUAW)/kBT^ + Nks 
7g = ^ • (60) 

In the limit of strong correlation one can replace (AUAW) with j{{AU)'^). Writing the resulting expression in terms of the 
excess (configurational) specific heat Cy gives 

7Cf? + NkB ,,,, 
7G = -y^ • (61) 

In the harmonic approximation, good for many simple liquids close to their melting point,^ Cy — {3/2)NkB (while it is 
generally larger in the supercooled liquid state). Thus the term Nks in the numerator is expected to be roughly a factor of ten 
smaller than the other term; we drop it and arrive at 

^-3^. (62) 

7 Cy 

This ratio is around one half in the harmonic approximation, otherwise larger. 



B. Energy-bond formulation of the strongly correlation property 

This section derives a general thermodynamic condition of the property of strong WU correlations, a condition which linearly 



constrains small variations in entropy, volume, temperature and pressure (Eq. (72i below). It is convenient to approach the 
problem from a general point of view. The ene rgy-bond formalism provides an abstract description of the interactions between 
a system and its surroundings An energy bond has an "effort" variable e{t) and a "flow" variable f{t), where 
e{t)f{t) is the free energy transferred into the system per unit time. The "displacement" q{t) is the time-integrated flow, i.e., 
q{t) = f{t). The energy-bond formalism is general, but we only discuss the linear case where it is most useful. Thus we consider 
a system that is slightly perturbed from equilibrium. It is assumed that the underlying microscopic dynamics is described by a 
stochastic equation, i.e., inertial forces are ignored. 

Linear-response theory is characterized by the fluctuation-dissipation (FD) theorem, which in the energy-bond formalism is 
given as follows. Consider a situation with n energy bonds and external control of the effort variables. If {fi{0)fj{t'))o is the 
equilibrium flow autocorrelation function, the average flow at time t is given by 

/^W=r^E/ {MO)m'))oej{t-t')dt'. (63) 
If the arbitrary additive constants of the displacements are chosen such that {qi)o = 0, the time-integrated version of this is 

-1 ^ /"OO 

I^W^-^T. {q^{0)fJ{t'))oe,{t-t')dt'. (64) 
^ j=i "'0 

If the flow variables are externally controlled, the FD theorem is 
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^^W^r^E/ (e.(0)e,(t')>o/,(i-0'^t'- (65) 

In most cases efforts are invariant under time reversal and flows change sign. The Onsager reciprocity relation is {fi{0)fj{t))o — 
{fj{0)fi{t))o (or {ei{0)ej{t))o = (e^ (O)ei(r))o, depending on which variables are externally controlled and which are free to 
fluctuate). From the FD theorem expressions for the frequency-dependent response functions are easily derived. Consider for 
instance the compliances Jij{uj), defined by Jij{uj) ~ 6qi{uj) /Sej{u!) for a periodic situation with linear perturbations around 
equilibrium, e{t) — Re[e(w) exp(ia;t)], etc. For these quantities the FD theorem implies 



M^) ^ {q^{0)f,{t'))oexp{-^iut')dt' . (66) 

f^B-^ Jo 

The case relevant to strongly correlating liquids is that of two energy bonds which are not independent, as we now proce ed to 
show. The two energy bonds are those of standard thermodynamics, reflecting the fundamental relation dE = TdS — pdyji§I12l 
The thermal energy bond with temperature variation as the effort variable and entropy variation as the displacement variable 
(ei{t) = ST{t), qi{t) = 5S{t)), and the mechanical energy bond with pressure variation as the effort variable and the negative 
volume variation as the displacement variable (e2{t) — Sp{t), q2{t) = —6V{t)). Usually the two standard thermodynamic 
energy bonds are independent, but we are here interested in the case when they are not. 

Treating the problem of two constrained energy bonds from a general perspective, we shall prove that the following four 
criteria are equivalent: 

1. The variables of the two energy bonds are Unearly constrained as follows 

aqi{t) +bq2it) = cei(t) + de2(i) . (67) 

2. The system's relaxing properties, i.e., its non-instantaneous responses, are describecPby a single variable e{t) as follows; 

qi{t) = Ji^ei(i) + Ji^e2(i) + 7ieW 

<?2(t) - J^,eiit) + J^e2{t) + ^2e{t). (68) 

In these equations the J°°'s are the compliances referring to the short-time, non-relaxing response (the high-frequency 
response). Note that = by the FD theorem. 



3. The relaxing parts of the three correlation functions entering into Eq. (66 1 are proportional. More precisely, the correlation 
functions obey 

(<7i(0)/i(t))o « (gi(0)/2(i))o « (g2(0)/i(t))o (X (g2(0)/2(i))o (MO), (69) 

and 

(9i(0)/i(i))o(g2(0)/2(i))o = ('Zi(0)/2(<))o('Z2(0)/i(i))o (MO). (70) 



Note that by differentiation Eq. ^ impHes for t ^ that (/i(0)/i(t))o oc (/i(0)/2(i))o oc (/2(0)/i(i))o oc 
(/2(0)/2(t))o. 

4. The dynamic Prigogine-Defay raticP^ A(a;) is unity at all frequencies (where double prime denotes the negative imaginary 
part): 



J'A{UJ)J^2H 



(71) 
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Proof that 1 2: By elimination of the variable e from Eq. (68 i, 2 implies 1. To prove the reverse implication, suppose 
that Eq. (67i applies and fix the dimensions such that the constants c and d are dimensionless. Define = (1 + c)/a, 
= J|f — ~l/b, and = {db + a)/6^. Introducing the variables ei = qi — Jif^i ~ J^S2 and £2 = Q2 ^ ^21^1 ^ J^^2, 
it follows that aei + be2 = 0. This means that we are in the situation described by Eq. (68 1 with a common relaxing variable to 
the two energy bonds, e{t) (x ei{t) cx e2(i), and symmetric short-time compliances, J{f — J|^. 

Proof that 2 ^ 3: In terms of functional derivatives with respect to the efforts at an earlier time {t' < t), since e{t) for small 
variations in the effort variables is linear in these, via the FD theorem time-reversal invariance implies that 5qi{t)/6e2{t') — 
6q2{t)/5ei{f). Thus Eq. (|68]l implies 5e{t)/5e2{t') cx 5e{t)/6ei{t'). From this Eqs. ^ and (|70|i now follow via the FD 
theorem. 

Proof that 3 4: According to the FD theorem the compliance matrix imaginary parts are given by Jij{oj) — 
(l/fcsT) /o°°('Zi(0)/j(i'))o sm{ujt')dt' . In conjunction with Eqs. (69 1 and (70 1 this implies that A{lj) = 1 at all frequencies. 

Proof that 4 2: We refer to the calculations of Ref. |48] which considered a system described by stochastic dynamics, i.e., 
with no inertial forces. Generalization of the arguments given there for the two standard thermodynamic energy bonds to the 
case of two arbitrary energy bonds proves the required implication. 

This completes the proof of the equivalence of points 1-4. - For the case where the two energy bonds are the fundamental 
thermodynamic bonds, the constraint Eq. (67 1 translates into (changing here the sign of b) 



aSS{t) +bdV{t) = cST{t) + dSp{t) . 



(72) 



How does this all relate to strong WU correlations in liquids? Via the equivalence of Eq. (72i to Eq. (681 and to unity 
dynamic Prigogine-Defay ratio (Eq. ([7T]i), the results derived in Refs. Ill2l3l implv that Eq. (72 1 describes a 100% correlating 
liquid subjected to small perturbations from equilibrium. Generally, for any strongly correlating liquid Eq. ( [72| is obeyed with 
good accuracy. Thus Eq. ( 72 1 gives the required thermodynamic formulation of the property of the hidden scale invariance 



characterizing strongly correlating liquids. 

Equation ^T2\ implies that for strongly correlating liquids the four thermodynamic variables, entropy, volume, temperature, 
and pressure cannot vary independently. Referring to Eq. ( 68 1, it is clear that for certain simultaneous changes of the four 



thermodynamic variables, the relaxing part is left unchanged; this suggests that for such changes the system is taken to a state 
where it is immediately in thermal equilibrium. This observation inspired the works leading to Paper IV where "isomorphs" 
are introduced. These are curves in the state diagram along which several quantities are invariant, and along which jumps from 
equilibrium at one state point take the system to a new state that is instantaneously in thermal equilibrium. 

Finally we would like to draw attention to an analogue of strongly correlating liquids. Consider a relaxing dielectric such 
as, e.g., a highly viscous dipolar liquid placed in a metal capacitor. This system's interaction with its surroundings may be 
described by two energy bonds: One energy bond is defined by the capacitor charge (electronic plus induced) and the voltage 
across the capacitor, the other energy bond is the induced dielectric charge at the capacitor surface and a fictive electric field only 
coupling to the liquid's dipoles. Because of Gauss' law these two energy bonds are not independent, but constrained by a linear 
displacement-field relation of the form Eq. ( [67] i. Thus from the energy-bond formalism point of view, a strongly correlating 
liquid is analogous to the standard measuring cell used for probing e{u!) of dipolar viscous liquids, with the strong virial / 
potential energy correlations reflecting one of Maxwell's four equations. 



VII. CONCLUDING REMARKS 



We have illuminated a number of features of strongly correlating liquids' hidden scale invariance. The linear term in the eIPL 
potential, which hides this approximate scale invariance, contributes little to the thermal fluct uations at fixed volume; this is why 
strongly correlating liquids inherit a number of IPL properties. As shown in previous paper^^HHls] jjjg bidden scale invariance 
has important experimental consequences, including that of density scaling. The general physical picture that we would like 
to suggest is that van der Waals liquids and most or all metallic liquids - because they are strongly correlating - are simpler than 
hydrogen-bonding liquids, ionic liquids, and covalently bonded liquids, which are not strongly correlating. 

Paper IV further investigates the consequences of a liquid being strongly correlating. This is done by defining "isomorphs" 
in the liquid's state diagram and showing that a number of properties to a good approximation are invariant along isomorphs. 
The isomorph definition does not refer to WU correlations. Only strongly correlating liquids have isomorphs, however; this is 
because the existence of isomorphs is a direct consequence of the hidden scale invariance of strongly correlating liquids. 
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APPENDIX A: CALCULATING THE ffcx, V CORRELATION COEFFICIENT IN THE NPT ENSEMBLE 



Here we provide the details of the calculation of the correlation coefficient between volume and excess enthalpy in the NpT 
ensemble. We apply Eq. (37 1 with {A, B} S {U, V}. First we need the pressure derivatives at constant temperature of ([/) and 
ly). Taking U first, we have (noting that for simple averages like ([/), it is not necessary to specify the ensemble because of 
equivalence of ensembles) 



d{U) 



dp 



d{V) 



dp 



d{U) 



dV 



Kt dV 



(Al) 
(A2) 



Note that V in the derivative is without averaging signs since there it is a parameter of the relevant ensemble (NVT). The volume 
derivative of (U) is calculated as follows. The excess (configurational) partition function Z{V, T) is the integral 



Z{V,T)^ J exp{^/3U{T,V))dT. 



(A3) 



Here F indexes points in configuration space and dT = d^^r/V^ . In the following we use the configuration space identity 
dU jdV = —WjV; constant temperature is implicit, as is the dependence of [/ on F and V , 



dV 



d_ 
dV 



j Uexp{-f3U)dT^ 



r 
dU_ 
dV 



dU 

Ui-p)—]expi-f3U)dT 



^ (^I^Uexp{-(3U)dr 



Z2 
dU 
W 

idU_ 
1 

V 



9Z 
dV 



BTT 

)-P{u—)-{u)z-' 



f3{U 



dV 
W 



{-f3)^exp{-l3U)dT 



p{AUA 



-m{ 



dU 



dU 



dV , 

{-{W) + I3{AUAW)) 



Thus we have (adding the subscript NVT to the fluctuation expression since this is ensemble-dependent) 



d{U) 



dp 



{W) (AC/AH/) NVT 



Kt 



(A4) 

(A5) 
(A6) 
(A7) 
(A8) 
(A9) 

(AlO) 



The pressure dependence of {V) is given by 



d{V) 



dp 



Kt' 



(All) 



To keep the notation simple, averaging signs are henceforth omitted from simple averages such as {V), {W), etc. We can write 



expressions for the variances of U and V in the NpT ensemble using Eq. (49 1 



keTKr f W {AUAW)nvt\ 



{{AV)' 



I NpT 



= 



V K^ 



V \Kt 
Kt 



kgTKT 



(A12) 
(A13) 
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We need also the covariance 

keTW {AUAW) 

- W + k:^ • ^^^^> 

Now we have all we need to construct the correlation coefficient in the NpT ensemble. The covariance between iJex 

and V is 

{AH,^AV)npt = {AUAV)npt + p{{AVf)NpT (A16) 
^ {AUAW)mvt keTW ^ VksT ^^^^^ 

{AUAW)nvt ^kBTiNksT) 
= + Kr ' ^^^'^ 

where we have used pV = NksT + W. The variance of ifex is more tedious: 

{{AH,^f)NpT = {{AUf)NpT + p'{{AVf)NpT + 2p{AUAV)NpT (A19) 

2 



/(Am^\ , knTKT ( W {AUAW)nvt \ 



+ + ^ {-keTW + {AUAW)nvt) (A20) 

/(Am^\ ^ / 2 2W{AUAW)nvt , {AUAW)},yr 

= {{^U) )nVT + TTIT- W r— + 



VKt \ ksT {keTf 

NVT \ 

~k^ ) 



Again using pV = NksT + W allows some simphcation: 



2 



((AH -((Am^\ ^^^^(w^ '^W{AUAW)mvt , {AUAW)^y^ 

((Aifex) )npT - ((AC/) ) +yj^[W ^ + (fc^r)2 



+ ^^(A^AJ)... ^ 2Ar(AC/AW^) ) (A22) 



+ W'^ + 2WNkBT + {NkeTf - 2W'^ - 2NkBTW 

.UA^ 
ksT 

-UAm^\ , {AUAW)%y^ N\kBTf ^NksT 
Now we can form the HeyV correlation coefficient, 

R = (AgexAV)^pT 24) 



{{AUAW) NVT + N{kBT)^)/KT 



VksT 
Kt 



V {{AU)^}nVT + kBTVKT + VKt + "iTT?^ (At/ AM/) 

{AUAW)nvt + NiksT)^ 

^KTVkBT{{AU)^) NVT + {AUAW)%yj. + N'^{kBT)'^ + 2N {AUAW) NVT (kBT)"^ 
a _ 1 



(A25) 

(A26) 
(A27) 
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where a = {AUAW)nvt + NikeTf and 6^ = KTVkBT{{AUf)NVT- 
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